I figured out the problem was that the expected mass in cylinders calculation was assuming PBC. I'm gonna have to actually compute the total mas in random cylinders. I'm gonna develop that here, as well as whatever I have to do for wp and cic.
In [1]:
import numpy as np
from glob import glob
from os import path
In [2]:
from matplotlib import pyplot as plt
%matplotlib inline
import seaborn as sns
sns.set()
import matplotlib.colors as colors
In [3]:
#from pearce.mocks.kittens import DarkSky
from pearce.mocks.kittens import TestBox
from halotools.empirical_models import Zheng07Cens, Zheng07Sats
from collections import OrderedDict
from time import time
import yaml
In [4]:
output_dir = './'# '/home/users/swmclau2/Git/pearce/bin/covmat/ds14_covmat/'
In [5]:
config_fname = 'xi_cosmo_trainer.yaml'
with open(path.join(output_dir, config_fname), 'r') as ymlfile:
cfg = yaml.load(ymlfile)
nd = float(cfg['HOD']['fixed_nd'] )
min_ptcl = int(cfg['HOD']['min_ptcl'])
r_bins = np.array(cfg['observation']['bins'] ).astype(float)
hod_param_ranges = cfg['HOD']['ordered_params']
logMmin_bounds = hod_param_ranges['logMmin']
del hod_param_ranges['logMmin']
In [6]:
def make_LHC(ordered_params, N, seed = None):
if seed is None:
seed = int(time())
np.random.seed(seed)
points = []
# by linspacing each parameter and shuffling, I ensure there is only one point in each row, in each dimension.
for plow, phigh in ordered_params.itervalues():
point = np.linspace(plow, phigh, num=N)
np.random.shuffle(point) # makes the cube random.
points.append(point)
return np.stack(points).T
def add_logMmin(hod_params, cat):
hod_params['logMmin'] = 13.0 #initial guess
#cat.populate(hod_params) #may be overkill, but will ensure params are written everywhere
def func(logMmin, hod_params):
hod_params.update({'logMmin':logMmin})
return (cat.calc_analytic_nd(hod_params, min_ptcl = min_ptcl) - nd)**2
res = minimize_scalar(func, bounds = logMmin_bounds, args = (hod_params,), options = {'maxiter':100}, method = 'Bounded')
# assuming this doens't fail
#print 'logMmin', res.x
hod_params['logMmin'] = res.x
In [7]:
def total_mass_enclosed_per_cylinder(centers, particles,
particle_masses, downsampling_factor, rp_bins, period,
num_threads=1, approx_cell1_size=None, approx_cell2_size=None):
# Perform bounds-checking and error-handling in private helper functions
#print period
args = (centers, particles, particle_masses, downsampling_factor,
rp_bins, period, num_threads)
result = _enclosed_mass_process_args(*args)
centers, particles, particle_masses, downsampling_factor, \
rp_bins, period, num_threads, PBCs = result
print 'A'
mean_particle_mass = np.mean(particle_masses)
normalized_particle_masses = particle_masses/mean_particle_mass
#print mean_particle_mass
#print normalized_particle_masses
# Calculate M_tot(< Rp) normalized with internal code units
total_mass_per_cylinder = weighted_npairs_per_object_xy(centers, particles,
normalized_particle_masses, rp_bins,
period=None, num_threads=num_threads, #try large finite PBCs
approx_cell1_size=approx_cell1_size,
approx_cell2_size=approx_cell2_size)
print 'B'
# Renormalize the particle masses and account for downsampling
total_mass_per_cylinder *= downsampling_factor*mean_particle_mass
return total_mass_per_cylinder
In [8]:
from halotools.mock_observables.surface_density.surface_density_helpers import rho_matter_comoving_in_halotools_units as rho_m_comoving
from halotools.mock_observables.surface_density.surface_density_helpers import annular_area_weighted_midpoints
from halotools.mock_observables.surface_density.surface_density_helpers import log_interpolation_with_inner_zero_masking as log_interp
from halotools.mock_observables import return_xyz_formatted_array
In [9]:
from halotools.mock_observables.surface_density.delta_sigma import _delta_sigma_precomputed_process_args
In [10]:
from halotools.mock_observables.surface_density.mass_in_cylinders import _enclosed_mass_process_args
from halotools.mock_observables.surface_density.weighted_npairs_per_object_xy import weighted_npairs_per_object_xy
In [11]:
def delta_sigma(galaxies, mass_enclosed_per_galaxy,
mass_enclosed_per_random, period,
rp_bins, cosmology):
# Perform bounds-checking and error-handling in private helper functions
args = (galaxies, mass_enclosed_per_galaxy, rp_bins, period)
result = _delta_sigma_precomputed_process_args(*args)
galaxies, mass_enclosed_per_galaxy, rp_bins, period, PBCs = result
total_mass_in_stack_of_cylinders = np.sum(mass_enclosed_per_galaxy, axis=0)
total_mass_in_stack_of_annuli = np.diff(total_mass_in_stack_of_cylinders)
mean_rho_comoving = rho_m_comoving(cosmology)
mean_sigma_comoving = mean_rho_comoving*float(period[2])
expected_mass_in_random_stack_of_cylinders = np.sum(mass_enclosed_per_random, axis = 0 )
expected_mass_in_random_stack_of_annuli = np.diff(expected_mass_in_random_stack_of_cylinders)
one_plus_mean_sigma_inside_rp = mean_sigma_comoving*(
total_mass_in_stack_of_cylinders/expected_mass_in_random_stack_of_cylinders)
one_plus_sigma = mean_sigma_comoving*(
total_mass_in_stack_of_annuli/expected_mass_in_random_stack_of_annuli)
rp_mids = annular_area_weighted_midpoints(rp_bins)
one_plus_mean_sigma_inside_rp_interp = log_interp(one_plus_mean_sigma_inside_rp,
rp_bins, rp_mids)
excess_surface_density = one_plus_mean_sigma_inside_rp_interp - one_plus_sigma
return excess_surface_density
In [12]:
def calc_ds(cat, rp_bins, randoms, tm_gal = None, tm_rand = None):
n_cores = 4
x_g, y_g, z_g = [cat.model.mock.galaxy_table[c] for c in ['x', 'y', 'z']]
pos_g = return_xyz_formatted_array(x_g, y_g, z_g, period=cat.Lbox)
print pos_g.shape
x_m, y_m, z_m = [cat.halocat.ptcl_table[c] for c in ['x', 'y', 'z']]
pos_m = return_xyz_formatted_array(x_m, y_m, z_m, period=cat.Lbox)
pos_r = randoms
print 'A'
if tm_gal is None:
tm_gal = total_mass_enclosed_per_cylinder(pos_g/cat.h, pos_m/cat.h,
cat.pmass/cat.h, 1./cat._downsample_factor, rp_bins, cat.Lbox/cat.h,
num_threads=n_cores)
print 'B'
if tm_rand is None:
tm_rand = total_mass_enclosed_per_cylinder(pos_r/cat.h, pos_m/cat.h,
cat.pmass/cat.h, 1./cat._downsample_factor, rp_bins, cat.Lbox/cat.h,
num_threads=n_cores)
print 'C'
return delta_sigma(pos_g / cat.h, tm_gal, tm_rand,
cat.Lbox / cat.h, rp_bins, cosmology=cat.cosmology)/(1e12), tm_gal, tm_rand
In [13]:
cat = TestBox(boxno = 0, realization = 0, system = 'sherlock')
cat.load(1.0, HOD='zheng07', particles = True, downsample_factor = 1e-2)
In [14]:
# TODO seed here for constant HODs
# TODO maybe just do 5, 10 may be overkill
N = 10
LHC = make_LHC(hod_param_ranges, N, 24)
hod_dicts = [dict(zip(hod_param_ranges.keys(), vals)) for vals in LHC]
In [15]:
cat.populate(hod_dicts[1])
In [16]:
rp_bins = np.logspace(-1.0, 1.6, 19)
In [17]:
N_rand = len(cat.model.mock.galaxy_table)
randoms = np.random.rand(int(N_rand), 3)*cat.Lbox
In [18]:
N_rand
Out[18]:
In [21]:
def _expected_mass_enclosed_in_random_stack_of_cylinders(num_total_cylinders,
Lbox, rp_bins, mean_rho_comoving):
cylinder_volumes = Lbox*np.pi*rp_bins**2
print cylinder_volumes
expected_mass_in_random_cylinder = mean_rho_comoving*cylinder_volumes
print expected_mass_in_random_cylinder
print num_total_cylinders
return expected_mass_in_random_cylinder*num_total_cylinders
In [22]:
from halotools.mock_observables.surface_density.surface_density_helpers import rho_matter_comoving_in_halotools_units as rho_m_comoving
In [23]:
mean_rho_comoving = rho_m_comoving(cat.cosmology)
In [24]:
em_rc = _expected_mass_enclosed_in_random_stack_of_cylinders(tm_rand.shape[0], cat.Lbox/cat.h,\
rp_bins, mean_rho_comoving)
In [25]:
em_rc/np.sum(tm_rand, axis = 0)
In [26]:
1/(cat.h**2)
Out[26]:
In [27]:
tm_rand.shape[0]
In [28]:
x_m, y_m, z_m = [cat.halocat.ptcl_table[c] for c in ['x', 'y', 'z']]
pos_m = return_xyz_formatted_array(x_m, y_m, z_m, period=cat.Lbox)
In [ ]:
randoms = np.random.rand(int(1e7), 3)*cat.Lbox
In [ ]:
tm_rand_big = total_mass_enclosed_per_cylinder(randoms/cat.h, pos_m/cat.h,
cat.pmass/cat.h, 1./cat._downsample_factor, rp_bins, cat.Lbox/cat.h,
num_threads=4)
In [ ]:
np.save('~/Git/pearce/bin/covmat/ds_14_covmat_v2/total_mass_randoms.npy', tm_rand_big)
In [ ]: